
# ==============================================================================
# DESCRIPTION ----
# ==============================================================================

# This code explores descriptive statistics from the cleaned amakudari.RData
# and outputs tables and figures. 

# ==============================================================================
# LIBRARIES AND IMPORT ----
# ==============================================================================

# Options
options(scipen = 100)

# Libraries
library(tidyverse)
library(igraph)
library(circlize)
library(viridis)
library(lubridate)
library(scales)
library(ggridges)

# Data import
amakudari <- read_csv("data/amakudata.csv")

# Line break function
breaks = function(x){gsub("([^ ]+ [^ ]+) ", "\\1\n", x)}

# ==============================================================================
# NUMBER OF HIRES BY FIRMS ----
# ==============================================================================

# Sum total hires by firm
hires_by_firm <- amakudari %>%
  mutate(large = ifelse(firm_dest_kana == "99", 1, 0)) %>%
  group_by(firm_dest_clean, large) %>%
  summarize(hires = n())

hires_by_firm_group <- hires_by_firm %>%
  group_by(hires) %>%
  summarize(firms = n())

# FIGURE A1: Plot distribution as empirical cumulative distribution function
total_hires_ecdf = 
  ggplot(hires_by_firm, aes(hires, group = large)) +
  stat_ecdf(geom = "smooth", aes(color = as.factor(large)), size = 0.75) +
  scale_color_manual(name = "Type of firm",
                     labels = c("Large", "All other"), 
                     values = c("seagreen3", "steelblue2")) +
  geom_vline(xintercept = 1, linetype = "dashed") +
  scale_x_continuous(limits = c(0, 15), breaks=seq(0, 15, 1)) +
  scale_y_continuous(limits = c(0.6, 1), breaks=seq(0.6, 1, 0.05),
                     labels = scales::percent_format(accuracy = 1)) +
  xlab("Number of hires per firm") + 
  ylab("% of firms with X number of hires or fewer") + 
  theme_bw() +
  theme(panel.grid.major = element_line(colour = "grey97")) +
  theme(panel.grid.minor = element_line(colour = "grey97")) +
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(axis.text=element_text(size = 10)) +
  theme(axis.text.x = element_text(size = 10)) +
  theme(legend.position = "bottom") +
  theme(legend.key.width = unit(3, "line"))

ggsave("figures/a1_hires_by_firms.pdf", height = 4, width = 7)

# ==============================================================================
# NUMBER OF AMAKUDARI BY MINISTRY ----
# ==============================================================================

#### Allow for adjustment for number of employees in each ministry ####
# Read in Japan Statistical Yearbook data
employees <- read_csv("data/jsy_agency_employees.csv") %>%
  filter(!is.na(ministry_short)) %>%
  select(-agency) %>%
  pivot_longer(!ministry_short, names_to = "year", values_to = "employees")

# Create count of hires by ministry and year
ministry_year <- amakudari %>% 
  mutate(year = as.character(year(date_ret))) %>%
  filter(ministry_short != "-99") %>%
  group_by(ministry_short, year) %>%
  summarise(Total = n()) %>%
  arrange(-Total)

ministry_year <- left_join(ministry_year, employees)

ministry_employees <- ministry_year %>% 
  filter(ministry_short != "-99", year != "2019", !is.na(year), 
         !is.na(employees), ministry_short != "MOD") %>% # Exclude MOD because of how employee numbers are recorded
  group_by(ministry_short) %>%
  summarise(hires = mean(Total), employees = mean(employees)) %>%
  mutate(percent = 100 * (hires/employees)) %>%
  arrange(-percent)

#### Adjusted chart of hires by ministry #### 
ministry <- amakudari %>% 
  filter(ministry_short != "-99") %>%
  group_by(ministry_short) %>%
  summarise(Total = n()) %>%
  arrange(-Total)

# FIGURE A3a: total amakudari appointments by former ministry
ggplot(ministry, aes(x = reorder(ministry_short, Total), y = Total)) +
  geom_bar(stat = "identity", width = 0.025, fill = "steelblue2") +
  geom_point(stat = "identity", shape = 21, size = 2, 
             fill = "steelblue2", color = "steelblue2") +
  geom_text(stat = 'identity', aes(label = Total), nudge_y = 100, 
            size = 2, color = "grey50") +
  scale_y_continuous(limits = c(0, 3000), breaks=seq(0, 3000, 500)) +
  xlab("") + 
  ylab("Count of hires by ministry") + 
  theme_classic() + 
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(axis.text=element_text(size = 10)) +
  theme(axis.text.x = element_text(size = 10)) + 
  coord_flip()

ggsave("figures/a3a_hires_by_ministry.pdf", height = 2.5, width = 6)

# FIGURE A3b: Mean amakudari appointments as a percentage of ministry employees
ministry_employees %>%
  filter(ministry_short != "MOD") %>%
  ggplot(aes(x = reorder(ministry_short, percent), y = percent)) +
  geom_bar(stat = "identity", width = 0.025, fill = "steelblue2") +
  geom_point(stat = "identity", shape = 21, size = 2, 
             fill = "steelblue2", color = "steelblue2") +
  geom_text(stat = 'identity', aes(label = round(hires, 0)), nudge_y = 0.1, 
            size = 2, color = "grey50") +
  scale_y_continuous(limits = c(0, 2.5), breaks=seq(0, 2.5, 0.5)) +
  xlab("Ministry") + 
  ylab("Percentage of ministry employees joining amakudari positions") + 
  theme_classic() + 
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(axis.text=element_text(size = 10)) +
  theme(axis.text.x = element_text(size = 10)) + 
  coord_flip()

ggsave("figures/a3b_hires_by_ministry_adjusted.pdf", height = 2.5, width = 6)

# ==============================================================================
# BAR CHART OF FIRM TYPES ----
# ==============================================================================

firm_type <- amakudari %>% 
  # Excluding NAs, including 99s
  filter(is.na(firm_type2_en) == FALSE,
         firm_type2_en != "99" ) %>% 
  mutate(firm_type2_en = ifelse(firm_type2_en == "na_public_corp", 
                                "Public Corporation", firm_type2_en), 
         firm_type2_en = ifelse(firm_type2_en == "na_government", 
                                "Government", firm_type2_en),
         firm_type2_en = ifelse(firm_type1_en == "Local government", # Change this in amakudata
                                "Government", firm_type2_en),
         firm_type4_en = ifelse(firm_type1_en == "Local government", 
                                "Government", firm_type4_en),
         dokuritsu = ifelse(str_detect(firm_dest, "独立行政法人"), 1, 0)
         )

# Change order of bars
positions <- c("Intermediary", "Government", "Public Corporation",
               "Non-stock", "Stock", "Public Interest" )

# Percentage public interest figure
firm_type %>%
  mutate(firm_type2_en = ifelse(str_detect(firm_type2_en, "Public"),
                                "Public & Public Interest", firm_type2_en)) %>%
  group_by(firm_type2_en) %>%
  summarise(n = n(), Percentage = (n() / nrow(firm_type)) * 100)

# FIGURE 1: Create plot
ggplot(firm_type, aes(firm_type2_en)) +
  geom_bar(aes(fill = firm_type4_en), stat = "count", width = 0.025) +
  geom_point(aes(fill = firm_type4_en, color = firm_type4_en), 
             stat = "count", shape = 21, size = 2) +
  geom_text(stat='count', aes(label=..count..), 
            nudge_y = 300, size = 2, color = "grey50") +
  xlab("Firm type") + 
  ylab("Count of hires by firm type") + 
  scale_x_discrete(breaks = positions,
                   limits = positions) +
  scale_y_continuous(limits = c(0, 6000), breaks=seq(0, 6000, 500)) +
  scale_fill_manual(values = c("firebrick3", "seagreen3", "steelblue2")) +
  scale_color_manual(values = c("firebrick3", "seagreen3", "steelblue2")) +
  theme_classic() + 
  theme(
    axis.title.y=element_blank(),
    axis.text=element_text(size = 9),
    axis.text.x = element_text(size = 9),
    legend.position = "bottom",
    legend.title=element_blank(),
    legend.key.size = unit(0.75, "line"),
    #legend.margin=margin(0,0,0,0),
    legend.box.margin=margin(-10, 0, 0, 0)
  ) +
  coord_flip()

ggsave("figures/1_hires_by_firm_type.pdf", height = 2, width = 6)

# Position level and firm type -------------------------------------------------
# Identify top positions
top_positions <- c("Vice-Minister", "Assistant Vice-Minister", "Comissioner",
                   "Director General", "Deputy Director General", "Director",
                   "Assistant Director", "Managing director")

# Create variables identifying top positions
firm_type <- firm_type %>%
  mutate(
    top_positions = ifelse(position_level %in% top_positions, 1, 0),
    vm = ifelse(position_level %in% c("Vice-Minister", "Assistant Vice-Minister"),
                1, 0)
  )

# FIGURE 2: Vice ministers
firm_type %>%
  group_by(vm) %>%
  mutate(total = n()) %>%  # calculate total count of each top_position
  mutate(firm_type2_en = ifelse(str_detect(firm_type2_en, "Public"),
                                "Public & Public Interest", firm_type2_en)) %>%
  group_by(vm, firm_type2_en) %>%
  summarise(n = n(), Percentage = (n() / first(total)) * 100) %>%
  select(-n) %>%
  ggplot(aes(fill=as.character(vm), y=Percentage, x = firm_type2_en)) + 
  geom_bar(position="dodge", stat="identity", width = 0.5) +
  scale_fill_manual(name = "Position level",
                    labels = c("All other positions", 
                               "Vice and assistant vice ministers"), 
                    values = c("steelblue2", "seagreen3")) +
  geom_text(position=position_dodge(width = 0.5), stat = 'identity', 
            aes(label = paste0(round(Percentage, 0), "%")), vjust = -2,
            size = 3, color = "grey50") +
  scale_y_continuous(limits = c(0, 60), breaks=seq(0, 60, 10)) +
  xlab("") +
  ylab("") +
  theme_classic() +
  theme(legend.position = "bottom")

ggsave("figures/2_vm_firm_dest.pdf", height = 4, width = 6)

# ==============================================================================
# FIRM TYPE BY MINISTRY ----
# ==============================================================================

# FIGURE A6a
firm_type %>%
  group_by(ministry_short) %>%
  mutate(total = n()) %>% 
  filter(ministry_short != -99) %>%
  group_by(ministry_short, firm_type2_en) %>%
  summarise(n = n(), Percentage = (n() / first(total)) * 100) %>%
  select(-n) %>%
  ggplot(aes(y=Percentage, x = ministry_short)) + 
  geom_bar(position="dodge", stat="identity", width = 0.001, 
           fill = "steelblue2", color = "steelblue2") +
  geom_point(stat = "identity", shape = 21, size = 2, 
             fill = "steelblue2", color = "steelblue2") +
  geom_text(position=position_dodge(width = 0.5), stat = 'identity', 
            aes(label = round(Percentage, 0)), vjust = -1,
            size = 3, color = "grey50") +
  scale_y_continuous(limits = c(0, 100), breaks=seq(0, 100, 10)) +
  facet_wrap(~firm_type2_en, nrow = 5) +
  xlab("") +
  theme_classic() +
  theme(legend.position = "bottom")

ggsave("figures/a6a_ministry_firm_type.pdf", height = 6, width = 10)

# # FIGURE A6b: exclude METI patent office workers
firm_type %>%
  filter(agency != "Patent Office") %>%
  group_by(ministry_short) %>%
  mutate(total = n()) %>% 
  filter(ministry_short != -99) %>%
  group_by(ministry_short, firm_type2_en) %>%
  summarise(n = n(), Percentage = (n() / first(total)) * 100) %>%
  select(-n) %>%
  ggplot(aes(y=Percentage, x = ministry_short)) + 
  geom_bar(position="dodge", stat="identity", width = 0.001, 
           fill = "steelblue2", color = "steelblue2") +
  geom_point(stat = "identity", shape = 21, size = 2, 
             fill = "steelblue2", color = "steelblue2") +
  geom_text(position=position_dodge(width = 0.5), stat = 'identity', 
            aes(label = round(Percentage, 0)), vjust = -1,
            size = 3, color = "grey50") +
  scale_y_continuous(limits = c(0, 100), breaks=seq(0, 100, 10)) +
  facet_wrap(~firm_type2_en, nrow = 5) +
  xlab("") +
  theme_classic() +
  theme(legend.position = "bottom")

ggsave("figures/a6b_ministry_firm_type_patent.pdf", height = 6, width = 10)

# ==============================================================================
# MINISTRY HIRES BY INDUSTRY ----
# ==============================================================================

# Create dataframe of top 10 retirement industries
top_industry <- amakudari %>%
  filter(
    industry != "99", 
    industry != "na_public_corp", 
    industry != "Public Interest",
    industry != "na_government",
    industry != "Intermediary",
    industry != "Others",
    firm_type1_en == "Stock company") %>%
  mutate(
    industry = case_when(
    industry == "Land Transportation" | industry == "Air Transportation" |
      industry == "Marine Transportation" ~ "Transportation",
    industry == "Securities & Commodity Futures" | 
      industry == "Other Financing Business" ~ "Finance",
    industry == "Nonclassifiable" & industry_detail == "Credit & Leasing" ~
      "Finance",
    TRUE ~ industry)) %>%
  group_by(industry) %>%
  summarise(total = n()) %>%
  arrange(-total) %>%
  slice(1:10) # Change to include more industries

# Create map of number of retirees from each ministry to each industry
ministry_industry_map <- amakudari %>%
  filter(
    ministry_short != "99", # Classified ministries only
    ministry_short != "-99") %>%
  mutate(
    industry = case_when(
    industry == "Land Transportation" | industry == "Air Transportation" |
      industry == "Marine Transportation" ~ "Transportation",
    industry == "Securities & Commodity Futures" | 
      industry == "Other Financing Business" ~ "Finance",
    industry == "Nonclassifiable" & industry_detail == "Credit & Leasing" ~
      "Finance",
      TRUE ~ industry)) %>%
  group_by(ministry_short, industry) %>%
  summarise(total_ministry = n()) %>%
  arrange(-total_ministry) %>%
  filter(!is.na(ministry_short)) %>%
  mutate(ministry = recode(ministry_short,
                           "MLIT" = "Land, Infrastructure, and Transport",
                           "MOF" = "Finance",
                           "MAFF" = "Agriculture, Forestry, and Fisheries",
                           "MIAC" = "Internal Affairs and Communications",
                           "MOJ" = "Justice",
                           "METI" = "Economy, Trade, and Industry",
                           "MOE" = "Environment",
                           "MHLW" = "Health, Labor, and Welfare",
                           "MOFA" = "Foreign Affairs",
                           "MOD" = "Defense",
                           "MEXT" = "Education, Sports, Science, and Tech"
  ))

# Merge together to create flow dataframe
top_industry_ministry <- 
  left_join(top_industry, ministry_industry_map,
            by = "industry") %>%
  select(-total) %>%
  select(ministry, industry, total_ministry) %>%
  group_by(industry) %>% 
  mutate(
    percent_total = total_ministry/sum(total_ministry),
    percent_total = percent(percent_total, accuracy = 1))

# FIGURE 3: Private sector hires by industry and ministry (2009 - 2019)
ggplot(top_industry_ministry, 
       aes(y = total_ministry, x = reorder(ministry, total_ministry),
           fill = ministry, color = ministry)) + 
  geom_bar(position="dodge", stat="identity", width = 0.005,
           fill = "steelblue2", color = "steelblue2") +
  geom_point(shape = 21, size = 1.5,
             fill = "steelblue2", color = "steelblue2") +
  geom_text(stat = 'identity', aes(label = percent_total), 
            nudge_y = 11, size = 2.75, color = "grey50") +
  facet_wrap(~ industry, nrow = 5) +
  scale_y_continuous(limits = c(0, 180), breaks=seq(0, 180, 25)) +
  theme_classic() +
  theme(legend.position="none") +
  xlab("") + 
  ylab("Total number of amakudari hires by industry") + 
  coord_flip()

ggsave("figures/3_industry_ministry_facet.pdf", height = 9, width = 9)

# ==============================================================================
# FIGURE A4: CHORD DIAGRAM PUBLIC INTEREST ----
# ==============================================================================

top_public <- amakudari %>%
  filter(
    firm_type2_en == "Public Interest", # Public corportations only
    ministry_short != "-99", # Classified ministries only
    firm_dest_en != "99", # Remove unclassified firms (or becomes largest total)
    !grepl("自営業", firm_dest)) %>% # Remove self employed
  group_by(firm_dest_en) %>%
  summarise(total = n()) %>%
  arrange(-total) %>%
  slice(1:10)

# Create map of number of retirees from each ministry to each firm
ministry_firm_map <- amakudari %>%
  filter(
    ministry_short != "99", # Classified ministries only
    ministry_short != "-99", # Classified ministries only
    firm_dest_en != "99",
    !grepl("自営業", firm_dest)) %>% # Remove self employed
  group_by(ministry_short, firm_dest_en) %>%
  summarise(total_ministry = n()) %>%
  arrange(-total_ministry) %>%
  filter(!is.na(ministry_short))

# Merge together
top_public_ministry = left_join(top_public, ministry_firm_map,
                                by = "firm_dest_en")

# Line break function
breaks = function(x){gsub("([^ ]+ [^ ]+) ", "\\1\n", x)}
top_public_ministry$firm_dest_en = breaks(top_public_ministry$firm_dest_en)

# Create chord dataframe
chord = top_public_ministry %>%
  select(ministry_short, firm_dest_en, total_ministry) %>%
  mutate(firm_dest_en = str_to_title(firm_dest_en))

# Create chord chart
chordDiagram(chord, 
             transparency = 0.25,  
             fontsize(10),
             big.gap = 10,
             small.gap = 2)
circos.clear()

unique(chord$firm_dest_en)

# Colors of chords
grid.col = c(MOFA = "black", MEXT = "grey15", MOE = "grey25",
             MOF = "gold", MOJ = "red", MIAC = "springgreen1", 
             MHLW = "cyan1", MLIT = "grey65", MAFF = "green4",
             METI = "steelblue4", 
             `Industrial Property\nCooperation Center` = "grey10",
             `Japan Forest\nFoundation` = "grey10",
             `Japan Meteorological\nBusiness Support\nCenter` = "grey10",
             `Japan International\nTrainee And\nSkilled Worker\nCooperation Organization` = "grey10",
             `Japan Correctional\nAssociation` = "grey10",
             `Airport Environmental\nImprovement Foundation` = "grey10",
             `Tokyo University` = "grey10",
             `Civil And\nJudicial Affairs\nAssociation` = "grey10",
             `Japan Hydrographic\nAssociation` = "grey10",
             `Japan Construction\nInformation Center` = "grey10"
)

# With horizontal labels
chord = chordDiagram(chord, 
                     grid.col = grid.col, 
                     transparency = 0.3,
                     big.gap = 10,
                     small.gap = 2,
                     annotationTrack = "grid",
                     preAllocateTracks = list(track.height = max(strwidth(unlist(dimnames(chord))))),
                     annotationTrackHeight = uh(1, "mm"))

# we go back to the first track and customize sector labels
circos.track(track.index = 1, panel.fun = function(x, y) {
  circos.text(CELL_META$xcenter, CELL_META$ylim[1], 
              CELL_META$sector.index, 
              cex = fontsize(6),
              facing = "clockwise", 
              niceFacing = TRUE, 
              adj = c(0, 0.5))
}, bg.border = NA) # here set bg.border to NA is important

# we go back to the first track and customize sector labels
circos.track(track.index = 1, panel.fun = function(x, y) {
  circos.axis(h = -0.07,
              major.at = c(0, 15, 30, 45, 60, 75, 90, 105, 120, 135),
              labels = TRUE,
              labels.cex = 0.4,
              direction = "inside",
              labels.facing = "clockwise",
              labels.pos.adjust = c(TRUE, TRUE),
              major.tick = TRUE,
              major.tick.length = 0.05,
              minor.ticks = 3,
              col = "grey85",
              labels.col = "white")
}, bg.border = NA) # here set bg.border to NA is important

dev.copy(pdf,'figures/a4_pub_interest_chord.pdf', width=10, height=10)
dev.off()
circos.clear()

# ==============================================================================
# FIGURE A5: CHORD DIAGRAM FOR-PROFIT FIRMS ----
# ==============================================================================
top_private <- amakudari %>%
  filter(
    firm_type3_en == "For-profit", # Public corportations only
    ministry_short != "99", # Classified ministries only
    ministry_short != "-99", # Classified ministries only
    firm_dest_en != "99",
    !grepl("自営業", firm_dest)) %>% # Remove self employed
  group_by(firm_dest_clean, firm_dest_en) %>%
  summarise(total = n()) %>%
  ungroup() %>%
  arrange(-total) %>%
  slice(1:10)

# Create map of number of retirees from each ministry to each firm
ministry_firm_map <- amakudari %>%
  filter(
    ministry_short != "99", # Classified ministries only
    ministry_short != "-99", # Classified ministries only
    firm_dest_en != "99",
    !grepl("自営業", firm_dest)) %>% # Remove self employed
  group_by(ministry_short, firm_dest_en) %>%
  summarise(total_ministry = n()) %>%
  arrange(-total_ministry) %>%
  filter(!is.na(ministry_short))

# Merge together
top_private_ministry = left_join(top_private, ministry_firm_map,
                                 by = "firm_dest_en")

# Line break function
breaks = function(x){gsub("([^ ]+ [^ ]+) ", "\\1\n", x)}
top_private_ministry$firm_dest_en = breaks(top_private_ministry$firm_dest_en)

# Create chord dataframe
chord = top_private_ministry %>%
  select(ministry_short, firm_dest_en, total_ministry) %>%
  mutate(firm_dest_en = str_to_title(firm_dest_en))

# Create chord chart
chordDiagram(chord, 
             transparency = 0.25,  
             fontsize(10),
             big.gap = 10,
             small.gap = 2)
circos.clear()

unique(chord$firm_dest_en)

# Colors of chords
grid.col = c(MOFA = "black", MEXT = "grey15", MOE = "grey25",
             MOF = "gold", MOJ = "red", MIAC = "springgreen1", 
             MHLW = "cyan1", MLIT = "grey65", MAFF = "green4",
             METI = "steelblue4", 
             `Mitsui Sumitomo\nInsurance` = "grey10",
             `Sompo Japan\nInsurance` = "grey10",
             `Tokio Marine\n& Nichido\nFire In` = "grey10",
             `Nippon Express` = "grey10",
             `Nec` = "grey10",
             `Fujitsu` = "grey10",
             `Pasona Group` = "grey10",
             `Mitsubishi Heavy\nIndustries` = "grey10",
             `Sumitomo` = "grey10",
             `All Nippon\nAirways` = "grey10"
)

# With horizontal labels
chord = chordDiagram(chord, 
                     grid.col = grid.col, 
                     transparency = 0.3,
                     big.gap = 10,
                     small.gap = 2,
                     annotationTrack = "grid",
                     preAllocateTracks = list(track.height = max(strwidth(unlist(dimnames(chord))))),
                     annotationTrackHeight = uh(1, "mm"))

# we go back to the first track and customize sector labels
circos.track(track.index = 1, panel.fun = function(x, y) {
  circos.text(CELL_META$xcenter, CELL_META$ylim[1], 
              CELL_META$sector.index, 
              cex = fontsize(6),
              facing = "clockwise", 
              niceFacing = TRUE, 
              adj = c(0, 0.5))
}, bg.border = NA) # here set bg.border to NA is important

# we go back to the first track and customize sector labels
circos.track(track.index = 1, panel.fun = function(x, y) {
  circos.axis(h = -0.07,
              major.at = c(0, 15, 30, 45, 60, 75, 90, 105, 120, 135),
              labels = TRUE,
              labels.cex = 0.4,
              direction = "inside",
              labels.facing = "clockwise",
              labels.pos.adjust = c(TRUE, TRUE),
              major.tick = TRUE,
              major.tick.length = 0.05,
              minor.ticks = 3,
              col = "grey85",
              labels.col = "white")
}, bg.border = NA) # here set bg.border to NA is important

dev.copy(pdf,'figures/a5_private_chord.pdf', width=10, height=10)
dev.off()
circos.clear()

# ==============================================================================
# FIGURE A7: AGE BY MINISTRY ----
# ==============================================================================
# Reorder factors for plotting
fct_order <- 
  c("MAFF", "METI", "MEXT", "MHLW", "MIAC", "MLIT", "MOD", "MOE", "MOF",
    "MOFA", "MOJ", "Other")

age_plot <- amakudari %>%
  mutate(
    ministry_short = ifelse(ministry_short == "-99", "Other", ministry_short),
    ministry_short = factor(ministry_short, levels = rev(fct_order))
  )

# Create plot
age_plot %>%
  ggplot(aes(x = age_ret, y = ministry_short, fill = ..x..)) +
  geom_density_ridges_gradient(scale = 4, rel_min_height = 0.01, color = "grey50") +
  geom_vline(xintercept = 60, color = "grey20") +
  scale_y_discrete(expand = c(0, 0)) +
  scale_x_continuous(expand = c(0, 0)) +
  coord_cartesian(clip = "off") +
  theme_ridges(center_axis_labels = TRUE) +
  scale_fill_viridis_c(option = "viridis", begin = 0.4, end = 0.9) +
  xlab("Age") +
  ylab("Ministry") +
  theme(
    panel.grid.major.x = element_blank(),
    panel.grid.minor.x = element_blank(),
    legend.position = "none"
    )

ggsave("figures/a7_age_ridge.pdf", height = 5, width = 7)
